Monocyte analysis V10

---

Analysis of monocyte data (SCZ vs CTR): most advanced model

Load data

First step QC

Plotting the first QC results:

1. Density plot

2. Covariate correlation

#Covariate correlation
MEs_lin(new.covs,new.covs)

3. PC-covariate correlation

#PC-covariate correlation (nomal p)
PCA_cov_cor_R(imputed.covs, filt.df)

Next step: QC on vst-normalized expression values.

1. Inter-sample distances

ComplexHeatmap::pheatmap(cor(filt.df,filt.df), color = heatmap.color.code)

2. PCA

#Check influence of covariates on data variance after transformation
pca <- plotPCA.custom(as.matrix(trnsf.df), intgroup=c("status", "batch", "sex", "smoking"), 
                      ntop = 50000, returnData=TRUE, metadata=imputed.covs, pc.1 = 1, pc.2 = 2)
PoV <- round(100 * attr(pca, "percentVar"))
PCAplot(pca, "batch", PoV.df=PoV, colors=c('orange','lightblue','lightgreen',"yellow"), pc.1 = 1, pc.2 = 2)
PCAplot(pca, "status", PoV.df=PoV, pc.1 = 1, pc.2 = 2)

3. Variance partitioning

#Check influence of covariates on data variance after transformation
plotVarPart(vp)

Data correction

1. PC-covariate correlation

print("The model used is:")

[1] “The model used is:”

dds.2@design

~Mean.Per.Base.Cov. + Mapped + Exonic.Rate + Genes.Detected + rRNA.rate + RIN + BMI + age_clusters + batch + status

#PC-covariate correlation (nomal p)
PCA_cov_cor_R(imputed.covs, batch.rem)

2. PCA plots

PCAplot(pca.cor, "batch", PoV.df=PoV.cor, colors=c('orange','lightblue','lightgreen',"yellow"), pc.1 = 1, pc.2 = 2)
PCAplot(pca.cor, "status", PoV.df=PoV.cor, pc.1 = 1, pc.2 = 2)

DESeq2 anaylsis

We then go into the differential gene-expression analysis:

1. Overview of the results

out of 24285 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 85, 0.35% LFC < 0 (down) : 111, 0.46% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 5) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

As well as the number of differentially expressed genes at lfc </> -/+0.5 at padj < 0.1:

[1] 43

Interactive results table

1. Results

#DT::renderDT(data.frame(res.complex), "OH1.5A",scrollY=1000)

createDT(data.frame(res), "Status", scrollY=1000)

Plotting results

2. Volcano plot:

EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'padj',
                pCutoff = 0.1,
                FCcutoff = 2,
                labSize = 6,
                xlim = c(-5,5),
                ylim = c(-1,5),
                legendPosition = 'right')

3. Plotting top genes (lFC < -1 / > 1) on a heatmap:

#Plotting DEGs
ComplexHeatmap::pheatmap(batch.rem[rownames(batch.rem) %in% upreg.genes,], scale= "row",
                         cluster_rows = T, cluster_cols = T, annotation_legend = T, show_colnames = F, show_rownames = T,
                         legend = T, treeheight_row = 0, color = heatmap.color.code,
                         annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)

ComplexHeatmap::pheatmap(batch.rem[rownames(batch.rem) %in% downreg.genes,], scale= "row",
                         cluster_rows = T, cluster_cols = T, annotation_legend = T, show_colnames = F, show_rownames = T,
                         legend = T, treeheight_row = 0, color = heatmap.color.code,
                         annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)

Enrichment testing results

1. List enrichments

#DT::renderDT(data.frame(res.complex), "OH1.5A",scrollY=1000)
paste("Reminder: we have", length(upreg.genes), "upregulated and", length(downreg.genes), "downregulated genes.")

[1] “Reminder: we have 19 upregulated and 10 downregulated genes.”

createDT(DEG.enrich.res, "Enrichment", scrollY=1000)

ComplexHeatmap::pheatmap(DEG.enrich, cluster_rows = F, cluster_cols = F, annotation_legend = F, show_colnames = T,
                         legend_breaks= c(0,-log(0.05), 5, 8, 10), 
                         legend = T, color = colorRampPalette(RColorBrewer::brewer.pal(6,"RdPu"))(30), 
                         legend_labels=c("0","2.995", "5", "8","-log(q) \n\n\n"))

2. GO enrichments

dotplot(GO.up, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free")

dotplot(GO.down, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free")

Variance causa testing

1. Gene list variance disparaties between groups

for (i in colnames(gene.panels_mrgd)){
  print(i)
  dat <- batch.rem[rownames(batch.rem) %in% na.omit(gene.panels_mrgd[,i]),]
  plots <- PCA_cov_cor_R(imputed.covs, dat)
  list.pca <- plotPCA.custom(as.matrix(dat), intgroup=c("status", "batch", "sex", "RIN"), 
                             ntop = 50000, returnData=TRUE, metadata=imputed.covs, pc.1=1, pc.2=2)
  list.pov <- round(100 * attr(list.pca, "percentVar"))
  print(PCAplot(list.pca, "status", Shape="batch", PoV.df=list.pov, pc.1=1, pc.2=2))}

[1] “Patir_microglia” [1] “Astrocytes” [1] “IFNy” [1] “NfKB” [1] “IFN.panel” [1] “LPS” [1] “IL4” [1] “GC.chronic” [1] “GC.acute.down” [1] “GC.acute.up”

#{r var.2, echo=TRUE, fig.height = 4, fig.width = 5.5, fig.align = "center"} #for (i in colnames(gene.panels_mrgd)){ # print(i) # dat <- batch.rem[rownames(batch.rem) %in% na.omit(gene.panels_mrgd[,i]),] # list.pca <- plotPCA.custom(as.matrix(dat), intgroup=c("status", "batch", "sex", "RIN"), # ntop = 50000, returnData=TRUE, metadata=imputed.covs, pc.1=1, pc.2=2) # list.pov <- round(100 * attr(list.pca, "percentVar")) # print(PCAplot(list.pca, "status", Shape="batch", PoV.df=list.pov, pc.1=1, pc.2=2))} #

Signature interrogation

1. Signature expression among most lFC genes (< -0.5 | > 0.5)

for (i in colnames(gene.panels_mrgd)){
  panel <- panel.list[[i]]
  DE <- DE.list[[i]]
  df <- panel[DE$log2FoldChange < -0.5 | DE$log2FoldChange > 0.5,]
  df <- df[,order(as.numeric(colnames(df)), method="radix", decreasing=F)]
  print(i)
  ComplexHeatmap::pheatmap(as.matrix(df), scale= "row",
                           cluster_rows = T, cluster_cols = F, annotation_legend = T, show_colnames = F, show_rownames = T,
                           legend = T, treeheight_row = 0, color = heatmap.color.code,
                           annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)
}

[1] “Patir_microglia” [1] “Astrocytes” [1] “IFNy” [1] “NfKB” [1] “IFN.panel” [1] “LPS” [1] “IL4” [1] “GC.chronic” [1] “GC.acute.down” [1] “GC.acute.up”

(1A) Compound median lFC for each signature

#median lFC to point directions
ComplexHeatmap::pheatmap(-compoundlFC, color = heatmap.color.code)

2. GEX-Status interaction

ComplexHeatmap::pheatmap(cbind(GEX.stat.cor[,1]), cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = T,
                         legend_breaks= c(0,0.01,0.05,0.1,0.1), 
                         legend = T, color = colorRampPalette(RColorBrewer::brewer.pal(6,"RdPu"))(30),
                         legend_labels=c("0","0.01","0.05","0.1","Rho \n\n\n"),
                         display_numbers = as.matrix(matrix_pvalue[,1]))

ComplexHeatmap::pheatmap(cbind(GEX.stat.cor[,3]), cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = T,
                         legend_breaks= c(0,1,2.5,2.5), 
                         legend = T, color = colorRampPalette(RColorBrewer::brewer.pal(6,"RdPu"))(30),
                         legend_labels=c("0","1","2.5","beta \n\n\n"),
                         display_numbers = as.matrix(matrix_pvalue[,2]))

(2B) As differential PC expression

scatters.all

3. DEG gene program- selected gene program interaction based on PCs

clean.compoundScatters

(3A) Gene program-Gene program interaction based on PCs

clean.compoundScatters2

(3B) Data summarized as correlation on heatmap

ComplexHeatmap::pheatmap(rcorr(as.matrix(cbind(DEGscatter.df,"Status"=scatter.df[,1])), type="spearman")$r,
                         color = heatmap.color.code, fontsize_row = 8, fontsize_col = 8)

4. Does status influence these interactions?

ComplexHeatmap::pheatmap(z, cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = T,
                         legend_breaks= c(-2,-1,0,1,2,2), 
                         legend = T, color = heatmap.color.code,
                         legend_labels=c("-2","-1","0","1","2","z-score \n\n\n"),
                         display_numbers = P)